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Abstract 

We consider a linear chain of quantum harmonic oscillators, in which the number of the individual 
^^ oscillators is given by an arbitrary number A^, and each oscillator is coupled at an arbitrary strength 

Q K to its nearest neighbors ("intra-coupling") , as well as the two end oscillators of the chain are 






O 



X 
S3 



coupled at an arbitrary strength Cy to two separate baths at arbitrarily different temperatures, 
respectively. We derive an exact closed expression for the steady-state heat current flowing from a 
hot bath through the chain to a cold bath, in the Drude-Ullersma damping model going beyond the 
Markovian damping. This allows us to explore the behavior of heat current relative to the intra- 



F^ coupling strength as a control parameter, especially in pursuit of the heat power amplification. Then 

Q it turns out that in the weak-coupling regime (k, Cy <^ 1), the heat current is small, as expected, and 

O 

'~~' almost independent of chain length N ^ hence violating Fourier's law of heat conduction, which is 

►^ consistent to the result in [IJ obtained within the rotating wave approximation for the intra-coupling 

QQ as well as the Born-Markov approximation for the chain-bath coupling. Beyond the weak-coupling 

T— I regime, on the other hand, we typically observe that with increase of the intra-coupling strength 

the heat current is gradually amplified, and reaches its maximum value at some specific coupling 
strength Kj^ "resonant" to a given chain-bath coupling strength. Also, the behavior of heat current 
^ versus chain length appears typically in such a way that the magnitude of current reaches its 

maximum with A^ = 1 and then gradually decreases with increase of the chain length, being in fact 
almost A^-independent in the range of A^ large enough. As a result, Fourier's law proves violated 
also in this regime. 
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I. INTRODUCTION 

The study of heat transport through small-scale quantum systems has recently attracted 
considerable interest due to an increasing demand for an understanding of the fundamental 
limit and efficiency of energy harvesting from a thermal machine at the quantum level [2]. 
One of the fundamental physical quantities considered in this subject is the heat current in 
the (non-equilibrium) steady state flowing from a hot bath through the quantum object of 
interest to a cold bath. 

The steady-state heat flux has been believed, for a long while, to obey Fourier's law of heat 
conduction stating that the heat flux is proportional to the gradient of temperature along 
its path, explicitly expressed as J = — Kp VT |3]; here the proportional constant Kp denotes 
the heat conductivity of the system in consideration, which is, typically for bulk materials, 
independent of the system size A^ and its shape, so giving rise to JT" oc 1/A^. In their seminal 
work, however, Rieder et al. discovered [4J that the steady-state heat flux through a one- 
dimensional classical harmonic chain is given hy J (x AT and so independent of the chain 
length (representing a novel form of energy flow), which accordingly deviates from Fourier's 
law. Since then, the validity (or not) of Fourier's law has come under scrutiny in various 
classical (e.g., [5]-tSj) and quantum systems (e.g., [Tl I91-I13)). Over the last few decades, in 
fact, it has turned out that Fourier's law may be violated in low-dimensional lattices whereas 
there is evidence that Fourier's law is still valid even for some one-dimensional classical 
and quantum systems. Therefore it remains an open question to rigorously determine the 
system-size dependence of the heat current. 

The rigorous analysis of heat transport through a small-scale quantum object has been 
carried out more recently. One of the interesting works is the study, given in [llj, of steady- 
state heat current through a disordered harmonic chain coupled to two baths at different 
temperatures, which was discussed mainly numerically, but giving rise to no clear conclusion 
regarding the system-size dependence of the heat current. Next, there was another inter- 
esting analytical treatment of this topic by Asadian et ai, given in jU [13], with a concrete 
conclusion that the heat current is independent of the system size accordingly violating 
Fourier's law of heat conduction. This was explicitly discussed in a harmonic chain coupled 
to baths as well as in a chain of two-level systems coupled to baths. In this treatment, 
they applied the Lindblad master equation formalism (within the Born-Markov approxima- 



tion) as well as the rotating wave approximation neglecting all energy non-conserving terms 
induced by the intra-coupling and so considering only the hopping of a single excitation 
between two nearest neighboring chain elements. As such, their analysis is restricted to 
the weak-coupling regime both in the chain-bath coupling and in the intra-coupling, as is 
typically the case for most studies. Accordingly, no sufficiently large energy flow can be an- 
ticipated to obtain. In addition, they demonstrated interestingly that Fourier's law can be 
recovered with chain length A^ — t- oo by adding, into the original Lindblad master equation, 
a superoperator representing the (phonon-induced) dephasing appearing in the condensed 
matter system. Needless to say, however, this additional dephasing Lindbladian cannot be 
derived from the original Hamiltonian describing the coupled chain plus baths under our 
current consideration. 

Therefore we are now demanded to study the steady-state heat flux beyond the weak- 
coupling regime, which has so far remained unexplored. By looking into this open problem, 
it is possible to examine behaviors of the heat flow relative to the coupling strengths as 
control parameters. This examination could stimulate the possibilities for increasing the ef- 
ficiency of energy harvesting by providing the amplified heat fiow followed by some additional 
novel quantum control of thermodynamic processes. In fact, the effects of dissipative en- 
vironments due to the system-bath coupling, which are normally negligible in macroscopic 
systems, become "detrimental" to low-dimensional quantum objects, and so the resultant 
noise is a major challenging factor to the control of, e.g., NEMS systems, as well-known 
[14| . Consequently, this subject is worthwhile to pay attention to, not only from the view- 
point of challenge in the quantum statistical physics but also from the viewpoint of quantum 
engineering. 

In the present paper, we consider a linear chain of quantum harmonic oscillators coupled 
at an arbitrary strength to two separate baths at different temperatures ("quantum Brow- 
nian harmonic chain"), in which each individual chain element is intra-coupled at another 
arbitrary strength to its nearest neighbors. We intend to provide an exact closed expression 
for the steady-state heat current through the harmonic chain as our central result [cf. Eqs. 



(41)-(41")]. The treatment of this physical quantity with rigor is mathematically manage- 
able due to the linear structure of our system. We approach this open problem by applying 
the quantum Langevin equation formalism to the Caldeira-Leggett type Hamiltonian. By 
doing this, we can go beyond the aforementioned weak-coupling approach. Our result may 



be straightforwardly generalized into a model of heat transport through a three-dimensional 
harmonic chain beyond the weak-coupling regime in case that heat diffusion perpendicular 
to the direction of the heat flux is neglected. 

The general layout of this paper is the following. In Sec. [TT] we briefly review and refine 
the general results regarding the quantum Brownian harmonic chain to be needed for our 



discussion and derive an exact closed expression of the bath correlation function. In Sec. Ill 



we 



we rigorously introduce a formal expression of the steady-state heat current. In Sec. |IV 
apply this formal expression to the simplest case of chain length A^ = 1 and derive an exact 
closed expression of the heat current. This result will be used as a basis for our discussion 
of the subsequent cases. In Sec. M we give an exact expression of the steady-state heat 



current for N = 2. In Sec. VI the same analysis will be carried out for an arbitrary chain 



length A^ > 3, giving rise to our central result. Finally we give the concluding remarks of 



this paper in Sec. VII 



II. EXACT EXPRESSION FOR THE BATH CORRELATION FUNCTION 

The linear chain of quantum Brownian oscillators under consideration is described by the 
model Hamiltonian of the Caldeira-Leggett type [15j 

H = Hs + Hb-^^sbi + Hbj^-sbN ) (1) 



where the isolated chain of A^ coupled linear oscillators, denoted by "system", 

N f p2 ./r \ N-1 

and two surrounding baths coupled to the first and the last oscillators of the chain are given 
by 

Hbi-sbi = 

respectively. Each of the two baths can split into the isolated bath and the system-bath 
coupling such as 

./=1 \ '^ / ,,—1 ,,—1 
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where the subscript fi = 1,N. Here the constant Kj is the (positive- valued) intra-couphng 
strength between two nearest neighboring oscillators of the chain with N > 2, and the set 
of constants {c^} denotes the (positive-valued) coupling strength between chain and each 
bath. Then the total system denoted by H is assumed initially in the separable state given 
by p(0) = Ps(0) ® Pis-^ ® Pfjj^- The local density matrix Ps(0) is an (arbitrary) initial state of 
the isolated chain Hg only, and the density matrix p^ = exp(— /3^i/(, )/^/3 is the canonical 
thermal equilibrium state of the isolated bath Hf, , where /3^ = 1/(A;bT^) and the partition 
function Z^ . Without any loss of generality, the two bath temperatures are assumed to 
meet Ti > T^. 

We apply the Heisenberg equation of motion to the Hamiltonian given in M, which can 
straightforwardly give rise to Pj{t) = MQj{t) and then the quantum Langevin equation |11] 

MQ,{t) + M / rfr7(t-r)A,,Q,(T) + MC.uQkit) = OW , (2) 

where for the sake of simplicity in form, the Einstein convention is applied in dealing with 
the subscripts j,k = 1,2, ■■ ■ , A^. Here the damping kernel and the shifted noise operator 
(representing a fluctuating force) are explicitly given by 

1 ""' C2 



0(t) = |-M7(t) Qm + ^bAt)\Si, + |-M7(t) g^(o) + a^W \ S^j , ([2 

respectively, where the fluctuating force of either isolated bath p, 

^bAt) = y^Cy{ x^,i.(0) cos(w,,t) + -^^ — sm{ujyt) \ . 



u=l 



Then it is easy to verify that the average value (^^(t))^^ = Tr{a^(t) p^^^} vanishes at any 
bath temperature /3^, as required. Next, the diagonal matrix Aj^ = 5jk {5ik + i^jvfc) connects 
the damping kernel to the two end oscillators directly coupled to the two separate baths. And 
the tridiagonal matrix of the isolated chain, Cjn = {f^^ + {i^n + /^n-i)/^} ^jn — {i^n ^j~i,n + 
i^n-iSj+i,n)/M = Cnji whcrc wc let Kq = Kjv = 0. It is also worthwhile to point out that 
this symmetric matrix C of real numbers is positive-definite, which can easily be shown by 
v*Cv > for any non-zero vector v of real numbers [16j, thus all its eigenvalues being 
positive- valued, as physically required in ([2]). 



For a physically realistic type of the damping kernel 7(t), we employ the form 7d(t) = 
7o cUji e"'^'* * of the well-known Drude-UUersma model, where a cut-off frequency Ud and a 
damping parameter 7^ fl7]. In this model, the bath correlation function in symmetrized 
form, defined as K^{t — t') := {{^b^(t) , ^b^(t')}+)i3^/2 in terms of the anti-commutator 
{ , }+, reduces to |18j 

K;f\t - f) = ^^^^ r du ^^ coth (^ cosMt - t')} . (3) 

After some algebraic manipulations, every single step of which is provided in detail in Ap- 
pendix [A, we can derive an exact expression of the correlation function K/j, {t — t'), given 
by 



where we introduce the effective frequency u^ = 27r kj^T^/h as well as the Lerch function 



$(2;, l,f) = f~^ 2-^1(1,"^; 1 + "i^; -2) [cf. (A6)]. Here the singularities of cot (vry) at y = 1, 2, •• • 
disappear due to their cancelation with those of the two Lerch functions at the same points; 
behaviors of K\i (t) versus time and temperature are plotted in Figs, hi and ^ respectively. 
From this closed expression, we can straightforwardly recover, in the limit of uj^ — ?■ 00, its 
classical counterpart given by K^^ {t — t') = ksT^^M 'jdi't — t') , being well-known. In addition, 
by letting the cut-off frequency Ud ^ 00 corresponding to the Ohmic (or Markovian) damp- 
ing 7o(t) = 27o5(t), the classical white-noise correlation function K^° (t — t') immediately 



follows, too. In fact, Eq. (3i) will play critical roles in deriving analytical expressions of the 



steady-state heat current in Sees. |IVHVI[ 

III. INTRODUCTION OF THE STEADY-STATE HEAT CURRENT 

We consider the average heat current flowing from a hot bath at temperature Ti through 
the coupled harmonic chain to a cold bath at T^ in the steady state p^''''\t), which is given 
by limi_j.oo p(t) in the Schrodinger picture. Accordingly, the steady-state energy expectation 
value (ifs)^(ss) of the harmonic chain is required to remain unchanged with the time. 



By substituting this into the Liouville equation dp^^^^dt = [H , p^""^^] / ih and then applying 
the cychc invariance of the trace, we can easily arrive at the expression 



Tr{ji'^)pW} = -Tr{ji-)pW}, 



(5) 



(N) 



where we have two energy current operators denoted by Ji^ = [Hg, Hb-^^sbi]/^^ and Jo 



?(iv) 



[Hs, Hi,^_aij^\/ih. We identify the left-hand side of (pi) as the steady-state input heat current 
jTJn and the right-hand side as the output heat current ~J'}ut, due to the fact that Ti > Tj^, 
we assume that J7in > 0. Then both heat currents J7in and —J'}ut, by construction, 
correspond to the input power from bath 1 into the harmonic chain and the output power 
from the chain to bath A^, respectively. 

Now let us find an explicit expression of the steady-state heat current J7in . We first sub- 
stitute (li)-(lo) and (2i) into (b| and then rewrite the operator Pi ® Xi^,y as {Pi , a;i,^}+/2, 
which will give rise to 



J. 



(N) 



^JA,-M7(0)Qi + I]c,Xi,J 



(6) 



From the Langevin equation given in Q with its damping term rewritten by integration by 
parts, we can also find that for a single Brownian oscillator (i.e., with chain length N = 1) 
coupled directly to both separate baths, 

5^c.x^,.(t) = ^4(t) + ^ {^2 + 27(0)} gi(t)+evW-^ {6iW+a„w} , (r) 



where /x = 1, 1' (i.e., 1' ^ A^), while for A^ > 2, 



m 



Next we substitute ([T]) into ^ and then into (IS]), as well as apply the cyclic invariance of 
the trace and p^''^ = \iTa.t^aoU{t) p{Qi)U\t) where U(t) = exp^—iHt/h). This allows us to 
switch from the Schrodinger picture to the Heisenberg picture. Then the steady-state heat 
current turns out to be 



J. 



(1) 



1 



lim Tr 



4M t^oo 



Pi(t),MQi(t) + MniQiit) + ^,,{t) - ^,^Xt)\^m 



(1) 



(8) 



for A^ = 1. In the same way, we can also find the corresponding expression of J'^J indepen- 
dently, shown to be identical to (8) but with exchange of ^bi{t) and C,b-^^,{t)- Similarly, Eqs. 

7 



(|5|, ^ and (7i) allow us to finally obtain the heat current 



J, 



(N) 



lim Tr 

2M t^oo 



A(t),MQi(t) + {Mnl + K,)Q,{t) - K,Q2{t)\ p(0) 



r{N) 



for A^ > 2, as well as its counterpart J'^^t , being identical to (81) but with substitution 



of Ki — )■ Kjv-i and, for all remaining subscripts, (1,2) -^ iN,N — 1). As shown, the key 
elements to the steady-state heat current are explicit expressions of {Qi{t),Q2{t)iPi{t)} in 
the limit of t — )■ 00. We will below restrict our discussion of these expressions, for the sake 
of simplicity, mainly to the case of fii = ■ ■ ■ = VL^^ =: Vt and /ti = ■ ■ ■ = Hm-i ='■ k- 

To derive an explicit form of each individual oscillator Qj{t), we directly apply the Laplace 
transform to the Langevin equation (2). Let its Laplace transform Qj{s) := C{Qj{t)}{s), 
then giving rise to C{Qj{t)}{s) = sQj{s) — Qj{0) and C{Qj{t)}{s) = s^ Qj{s) — sQj{0) — 
Qj{0) [19j. Then we can easily obtain 

1 



where the Laplace-transformed fluctuating force 



M 



^6i(s) 5^1+^62(5) (5 



JN 



CbAs) 






s2 + u;2 



5X^,1.(0) + 



P/.,.(0) 



m. 



(9) 



&) 



and the symmetric tridiagonal matrix Bjk{s) = s^ 6jk + ^7(5) Aj^ + Cjk expressed in terms 
of the Laplace-transformed damping kernel 

.2 1 



7(s) 



s 
M 






rriu ul s"^ + ujI' 



#) 



In the Drude-Ullersma model, the damping kernel 7(5) — )■ 7d(s) = 'j^Ud/is + Ud) [15j . 
Therefore, the central task to be undertaken is the determination of an explicit form of the 
inverse matrix B^^{s) =: A{s), which will be performed below for individual chain lengths 

N. 



IV. STEADY-STATE HEAT CURRENT FOR THE CASE OF iV = 1 



We begin with the simplest case of A^ = 1, in which a single oscillator is coupled directly to 
two separate baths at different temperatures. As well-known, the matrix A{s) then reduces 
to Mx{s), where 

X{s) = jj ,, , J, ^^,, , (10) 



M s2 + fi2 + S7(s) 



8 



corresponding to the dynamic susceptibility in the frequency domain, given by x{^) ^ x{s) 



with s —7- —iu + O"*" |15j. In the Drude-UUersma model, Eq. (10) reduces to 

is + cod)/M 



Xd[s) 



hi(s) 



(lOi 



where hi{s) = s^ + Ud s'^ + {Q"^ + '~fo^d) s + fi^ w^, with all its coefficients being positive- 
valued. Accordingly, this cubic polynomial can be factorized as {s + Zo){s + Zi){s + Z2), 
where Iie{zo),Iie{zi),Iie{z2) > 0, through the symmetric relations 



Zo + Zi + Z2 = Ud , fi^ + 7o ^d = Zq {Zi + Z2) + Z1Z2 , fi^ Ud = Zo Zi Z2 : 

[zq,Zi,Z2) can equivalently be rewritten as (wo,^;o?7); where [20] 



(11) 



n' 



(wo)' 



Zq 



UJd 



zo + 7 



7c = 7 



zq {zo + 7) + (wo)' 



then these lead to Zi = 7/2 + iwi and Z2 = 7/2 — iwi with Wi = a/(wo)^ — (7/2)^. The 
parameters [zq, Zi, Z2) will be useful for a compact expression of the steady-state heat current 



r{l) 



J7J,; [cf. (19)]. In fact, these can be explicitly expressed in terms of (i7, Wf^, 7o) by the cubic 
formula, as well-known |21] . 

Now let us find an explicit expression of the single oscillator Qi (t) in the limit of t — )■ 00 by 
considering the equation of its Laplace transform Qi{s) given in (9). This can be efficiently 
carried out with the aid of the final value theorem of the Laplace transform |22]; it reads 
as limt_>oo /(t) = lims_>o S-F(s), where F{s) = C{f{t)}{s), upon condition that all poles of 
F{s), except s = 0, have negative real parts, i.e., if sF{s) is analytic on the imaginary axis 



and in the right half-plane. By noting from (lOi) the fact that this condition is met by Xd{s) 
and so lim^^o sXdi-s) = linis^o s'^ Xdis) = 0, we can easily obtain from ^ the expression 



t—>oo 



t— >oo 



limQi(t) = lim / dt'xd{t-t') {di(0 + 4,,(0 



(12) 



where the response function [23j 
Xdit) = C-'{xdis)} -- 



-20* 



'Zl t 



-22 i 



M 



{Zq - Zi) {Zi - Z2) {Z2 - Zo) 



(12i) 



For the sake of comparison with (12), it is also worthwhile to mention that both fluctuating 



forces C,bi{s) and C,b ,(s), as shown in (9i), do not meet the prerequisite for applying the final 



value theorem, though, and so it turns out that 



lim Qi{t) ^ \imsxd{s) \^bi{s) + ^b,,{s) 



Similarly, it appears from Pi{t) = MQi{t) and Xd(0) = that 

lim A(t) = M lim / dt' !^xd{t-t')] {4^(t')+ev(0| • (13) 

Now we are ready to derive an explicit expression of the steady-state current. Substituting 



(12) and (13) into (pj), we first obtain the steady-state expectation value 

(s=) 



Piit),Qiit) 



2M lim / dr I dr' Xdit -t)\ ^Y^it - r ) > x 



t— ^oo 



dt' 



[k^^\t-t') +Kf{T-T')] 



(14) 



Plugging subsequently into this expression both (12i) and (3i) with (A6), and then using 

m 










a" + a' + 2a 







ia" + a') ia' + a) (a + a") ' 



(15) 



where a = Ud oi nuju with n = 0, 1, 2, ■ ■ ■ , we can explicitly evaluate the double integral in 



(14), which turns out to vanish due to the symmetric structure of Xd{t) in (^q, Zi,Z2). Next, 
taking into account the fact that Qi{t) = Pi(t)/M and dxd{0)/dt = 1/M, we can also find 
that ({Pi(t) , Qi(t)}+)^""'' = 0. Then the formal expression of the heat current given in (8) 
is simplified as 



J, 



(1) 



AM 



(ss) 



(16) 



Pi{t) , ^bAt) - ^bAt) ^ 

Along the same line, we can also obtain the expression for J'}J, identical to (16) but with 



exchange of ^bi{t) and ^b ,i'^)i which immediately verifies that J'}J = — jTji indeed. In the 
equilibrium state given by /3i = /3i/, the heat current vanishes, as expected. 

After some algebraic manipulations, every single step of which is provided in detail in 



Appendix [Bj Eq. (16) finally reduces to the exact expression 

Zj ■ {^d — Zj) 






1) h'y{zoUJd + {woy} 



(1) 



Att 



j=0 



7- — : ) (, ~ ^ {^^i(^j-) - '^%(^^-) 

[Zj - Zj+i) [Zj - Zj+2) ^ 



(17) 



in terms of the parameters (wq, Zq, 7) given in (111), where j = j (mod 3). Here 



'^fiiZj) 



-vrcot(i^)+^(-/^)-^(^) ^(^)-V^ 



I3hzj 
2tt 



UJd + Z^ 



Ud — Zj 



(18) 



10 



where the digamina function ip{y) = d \nT{y)/dy |25| . With the help of ([TT|-([TTi) and 



(BlO)-(Bll) as well as the relation given by —n cot(y) +ip{—y/-K) = ^lj(y/n) +-K/y, Eq. (17) 



can be rewritten as a compact expression 



j; 



(1) 



^di' 



d lo 



2ti 



hi{uJd) 



+ E 



i^j ^ ^i+i)(^i ~ ^i+2)(^j + ^d) 



=0 ^"i 



(19) 



where ipp{y) '■= tp{Phy/27i), and g^ := —{ipfilud) + 27r /Phud}- Now we can also consider the 
semiclassical behavior of the heat current by expanding the diganima functions such that in 
the limit oi f3h -^ 0, 



J: 



(1) 



J. 



(1) 



/i'/3i/3i' 



(2fi2 + 7„a;d) + 



;i4/3i/3i.(/32 + /3i/3i, + /32) 



24 ' "" ' '"""' ' 25 -32 -5 

{2fi^ + 7.a;,(3fi2 + ^^cu, + c^2)} + O {{pKf)\ , 



X 



(20) 



expressed in terms of the original input parameters (fi, w^, 7^) only, the derivation of which 
is provided in Appendix |B|) . Here the leading term 



J. 



(1) 



7oW^ 



2fi2 + ^^^, + 2a;2 y^^ ^^, 



1 



&) 



(1) 



corresponds to the classical counterpart to J^^ . In the Ohmic limit ujd -^ 00, this classical 
value reduces to the well-known expression given by J^:^ = (7o/2) ■ (l//3i — l//3i')- Here we 
see that the heat current J'}^ is not directly proportional to the bath-temperature difference, 
especially in the low-temperature regime. As a result, it turns out that Fourier's law of heat 
conduction is not valid even for the case of A^ = 1. The behaviors of jTii versus the "hot- 
bath" temperature Ti are plotted in Figs. |3] and |4| where the "cold-bath" temperature Ti' 
is imposed in the low-temperature and the high-temperature regime, respectively; in the 
weak-coupling limit imposed by 7^ ^ Q, the low-magnitude heat current is observed indeed. 



V. STEADY-STATE HEAT CURRENT FOR THE CASE OF N 



Next we consider the case of A^ = 2. To efficiently proceed with the determination of an 
explicit form of the inverse matrix B~^{s), we ffist diagonalize the tridiagonal matrix B{s). 
To do so, we introduce the normal coordinates {Qj{s)} of the isolated chain Hs given in 

jk 



(li), which satisfy Qj{s) = OjkQk{s) and {0'CO)jk = ^l^jk with O* O = 00* = 1^ [11]. 



Eq. ([9]) is then rewritten as 



%fc(s)Qfc(s) = sQ.iO) + Q,(0) + —\^{s)Oi,+^is)0^,f , 



(21) 



11 



where the matrix 



(6*i3o).^(s) = {s^ + nl) 5,k + s^is) Uk - 5^0„,0„fcj 



m) 



being, in fact, of diagonal form for N = 2. For the case of f2i = ^2 =• ^ to be considered 
here, it easily turns out that ^1=^7 and ^2 = (fi^ + 2k/ M)^^"^ as well as 



O 



which is a simple constant matrix; for an isolated chain with A^ > 3, in comparison, it is 
straightforward to verify that the matrix elements Ojk are not mere constants but functions of 
s [26]. Then the inverse matrix !B~^(s) appears as a diagonal form with (53^^)ii(s) = Mxi{s) 




and {^~^)22{s) = Mx2{s), where Xni^) with yU = 1,2 correspond to x{s) given in (10) used 
for A^ = 1, with substitution of i7 — )■ 0^, respectively. 

We are now ready to apply to each of these two diagonal elements the same technique as 



for A^ = 1. Then it turns out, with the help of (21 3), that 

2 



limgi(t) 



Iimg2(t) 

i— >oo Z 



1 /"* r - 

- lim drV] Xi.i't - r) u"bi 
2}}^ /''T.X,it-r){i 

-'^ u=l 



-(-ira.(r) 
6.(r)-(-ira.(r)L 



(22a) 
(22b) 



r(2) 



To explicitly evaluate the formal expression of the heat current J7i„ in (|8pi|), we first focus 
on the steady-state expectation value 



{A(t),gi(t) 



(ss) 



— lim / dT dT'-{xiit-r)x2it-r')}x 
2 t^^Jo Jo ot 



[k^^\ 



T — T 



K. 



W, 



T — T 



(23) 



By means of the same technique as for the case of A^ = 1, we can straightforwardly show 
that this vanishes indeed. In the same way, it also turns out that {{Pi{t) , gi(t)}+)*-''°'* = 0. 
Similarly, the last needed expectation value 

(bb) 



{A(t),g2(t)} 



- lim / dr / dr' [x,{t - r) {dt X2{t - r')} 

2 i^°°Jo Jo 



(24) 



{dtXi{t-r)}x2{t-r')] KT\r-r') - K'^\r - r' 



id), 



.{d). 



7^0 
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can be evaluated in closed form (cf. Appendix O. Substituting this form into (81), we can 
immediately arrive at an exact expression of the steady-state heat current 

fss) *- 9 



J: 



(2) 



K 



Plit),Q2it) 



2 ^^{ul-z])-z,-{T,^{z,)-T,^{z,)} 



(25) 



E 

i=o 



{Zj - Zj+i) ■ {Zj - Zj+2) ■ hi2{Zj) 



Zj -^ z' ; hi2{zj) -^ huizj] 



where the parameters Zj^s are identical to {zq, Zi, Z2) given in (11) but with substitution of 



hui-z'j) = 0. As given in (19) 



Q — !■ Cli, and so are z'-^s with Q -^ 0,2- The function hii{z',-) := hi{z'-) with Q — )■ Cli, and 
huizj) := hi{zj) with Q — > ^2; by construction, hii{—Zj 
for A^ = 1, Eq. (|25|) can finally be rewritten as 

1 1 






(2) 



'M2 



(I2) 



/3i /32 



+ 



27rM 



X 



(26) 






zr{^0iizj)-ij^^{zj)} 

{zj - Zj+i) ■ {zj - Zj+2) ■ huizj) 



Zj -^ z\\ hi2izj) -^ hu{z'-) 



Here (I2) := ujj ■ /k/(2 A^), where /« = fi^ + k/M + uj and 

K \3 



2 4/'*^,o2 2 4 



m) + 



Now we apply to the exact expression given in (26) the same technique as provided for 



(20), in order to study its semiclassical behavior. Then it turns out that 

(2) \, h' /3i (32 V, h^ (3, (32 {(31 + (3,(32 + (3l) fi^ u^ 



J: 



(2) 



•J cl 



o{mf] 



, (27) 



12 /. 23.32.5 /, 

where f ^ := (2 VL^ + 7o i^d + ^"d) ' 1^/^ + VL^ + Vt^ ujd (7o + ^d)- Here, the leading term is given 
by the classical counterpart 



J. 



(2) ^ ftTX 
M2 



(12) 



1 1 



|7k) 



^(2) ;. 



The behaviors of J7in in the low-temperature and the high-temperature regime are plotted 
in Figs. [5] and [6} respectively; as demonstrated, they are consistent with the behaviors of 
J7il . Here the weak-coupling limit is imposed by k/M, 72 <^ ^^2. Along the same line as 
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(23)-(24), we can also obtain that 

A(t),Q2(t)}^^ = -^{A(t),gi(t) 
A(t),4w}^\ = -/{A(t),4(t) 
p2{t),Q,it)]\ = -/{A(t),g2(t) 

From this, it follows that jj^^ = - Jii^^ indeed. We see from (|27l-(|27f]) that J,^^^ -^ for 



(ss) 



(ss) 



(ss) 



(28a) 
(28b) 
(28c) 



K — )■ 0, and so there can be no sufficiently high output power in the weak-coupling regime, as 

(2) 

expected. Finally we remark for a later purpose that the steady-state heat current J7in was 
rigorously treated based on the two uncoupled normal modes, each of which was thoroughly 
studied for J7il already in Sec. 



IV 



VI. STEADY-STATE HEAT CURRENT FOR THE CASE OF iV > 3 



We first need to point out that the matrix 23 (s) given in (21a,) is not of diagonal form 

for A^ > 3 and neither is its inverse. Therefore, the technique provided for A^ = 2 cannot 

be applied any longer. Instead, we adopt a different approach to the determination of an 

explicit form of B^^{s) = A{s), developed in [271128]; given an A^ x A^ symmetric tridiagonal 

matrix 

/acO 0---00\ 

c b c ••■ 

c b c ■■■ 

0c 6 ■■■ 



B{s) 



■■■ c b c 
\000--- c a J 



(29) 



where 



.2 , „.-,„^ , o2 , '^ 5 := s2^fi2^2«: 



a ■= s^ + s7(s) + n^ + 



M 



M 



c :- 



K 

M 



(29i 



Let 6 
form, 



s'~f{s) + c, and so a = b + 6. Then its inverse is explicitly given by a symmetric 



Ajk{s) 



1 Numjfc(s;A^) 



c sin0 Den(s; A^) 



(30) 
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for j < k, where both numerator and denominator are 



Numjfc(s; A^) = [c {sin j0} - 6 {sin(j - 1)0}] x 

[S {sm{N - A:)0} - c {sin(A^ -k + 1)0}] 
Den(s; A^) = c^ {sin(A^ + 1)0} -2c6 {sin N(j)} + 5^ {sin(A^ - 1)0} 



respectively. Here, 



sm( 



— (b'-4c'f , COS. 
2c ^ 



b_ 



(30 



For s G M, the functions, sin0 and cos0 are rewritten as ^sinhr and coshr, respectively, in 
terms of a real number r = —icf) > 0. 

Next let us express the matrix elements Ajk{s) explicitly in terms of s. By using sin(n + 
1)0 = 2 (cos 0) (sin 720) — sin(ra — 1)0 and 

(n — Z/) TT 



smnc 



u=0 



n 



(cos d)Y (sm.d)Y '^ sin 



(31) 



as well as sin(?27r/2) = i {(—«)" — 'i"}/2 [25j, we can rewrite Eqs. (30i) and (30o) as 

(_2c)'=-J+i-^ 



Numjfe(s;Ar) 
Den(s; A^) 
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i{-2cY-'' 



- {F,(s) + 2,5F,_i(s)}{F^_,+i(s) + 2<5F^_fc(s)} (32a) 
{{a + 5)F,{s)-2{^-5^)F,_,{s)] , (32b) 



respectively. Here -F„(s) = G'^{s) — H"'{s), where 

G{s) := 6+(62_4c2)V2 ^ H{s) := b-ib^-AcY^\ 



(323) 



From this, we see that Fi{s) = Aic sin0, and G = —2ce ^'^ and H = — 2ce*'^, as well as 



F„(s) = -2i (-2c)" sin(n0) , F„+i(s)/F„(s) = -2c{sin(ri + l)0}/{sinri(; 



IH) 



Substituting into (30) the expressions given in (32a) and (32b) as well as 7(5) — )■ 7d(s) 
7o(X'rf/(s + Ud), we can explicitly obtain 
(S + UJd) (s^ 



Aii{s) 

Al2{s) = -C' 



^^-^ + cuds'^-' + 



(s + w,) (52^^-3 + u;rfs2^-4 + 



AoNis) 



-c 



N-2 



/i]v(s) 



Ain{s) 



s + Ud) ■hi{s)/h^{s) 



{s + LOdY/h^{s), (33) 



where the cubic polynomial hi{s) appears from hi{s) given in (lOi) with Vl — > (fi^ + zt/M)^/^. 
Here we have the (2A^ + 2)th-degree polynomial, hi^{s) = den(s; A^) ■ (s + Ud)^ = a2iv+2 • 
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g2N+2 j^2 ujd s^^"*""^ + ■ ■ ■ , where the leading coefficient a2N+2 = 1, and the factor den(s; N) = 
—4i ■ Den(s; A^)/{(— c)^^^ Fi(s)}; den(s; A^) is not a polynomial, due to the fractional form 
of the damping kernel 7d(s), and its completely explicit expression is provided in Appendix 



Dj With the aid of (32^) and (32i), it can also be shown that hf^{—UJd) 7^ 0. 

Owing to these explicit expressions of Ajk{s)., we are now in position to straightforwardly 
proceed to obtain Ajk{t) = C^^{Ajk{s)}{t) in time domain. By applying the initial value 
theorem given by limi^o/(^) = linis^oo sF{s) [22j, we can easily find that Ajk{0) = and, 
e.g., Aii(O) = lims_s.oo s {s^ii(s) — ylii(O)} = 1 as well as ^12(0) = 0, etc. It can also be 
verified that hj^{s) = Y[j=o ('^ + ^j) nieets the condition of 'Re(zj) > indeed (cf. Appendix 



M, as the cubic polynomial hi{s) in (lOi) did. Due to this fact, we are also allowed to apply 



the final value theorem, then giving rise to limf^oo Ajk{t) = limt_>oo Ajk{t) = to be needed 
below. 

Now we explicitly consider the formal expression of the steady-state heat current given 
in 



1) based on the above result. To do so, we begin with 

lim Qi(t) = j- lim / dr \Au{t - r) ■ 4,(r) + Ai^{t - r) ■ ^^(r)) 



lim QaW = T7 lim / dr \ A^iit - r) ■ ^,,(r) + A^^it - r) ■ ^,^(r) 
lim A(t) = lim [ dTdt\ A^^{t - r) ■ ^^Ar) + Ai^(t - r) ■ ^.^(t) 



limQi(t) = T^lim 



drdt Anit - t) ■ ai(r) + A^4t - r) ■ ^^^(r) + 



— \un\An{0)-^M + A,,{0) ■ ^,^{t) 



(34a) 
(34b) 
(34c) 

(34d) 



which can be found from the equation of Laplace-transform in ([9]). Substituting (34b) and 



(34c) into (81), we acquire the ffist expectation value 

t 



lim Tr 

i— >oo 



{A(t),g2(t)} ■p(o) 



lim 







dr / dr' [{dtAuit-T)}-A2iit-T')x 







<d) 



Krir - r') + {dtA^^it - r)} ■ A^^t - r') ■ K'^'ir - r') 



(d), 



^0. 



(35) 



This integral can be evaluated explicitly in the same way as in (24) valid for a chain with 



N = 2 only [cf. (38) and (40)]; the detail of this evaluating process is provided in Appendix 



16 



[E] , where we also verify the equahty 



hm f dr [ dr' {dt Au{t - t)} ■ A2i{t - t') ■ e-''^^-^'^ 
■^°° Jo Jo 

- hm [ dr [ dr' {dtAi^{t - r)} • A2^(t - r') ■ e""!"""'! . 
*-^°° Jo io 



(36) 



Here the constant a equals the cut-off frequency Ud or the effective frequency given by nui 
or noj,^, where n = 0, 1, 2, ■ ■ ■ . Next we consider the second steady-state expectation value, 



linii_i.oo Tr[{Pi(t), Qi(t)}+ ■ /5(0)] by applying the same technique to (34a) and (34c); with 



the help of ( 36 ) , this will straightforwardly give rise to 

lim [ dr [ dr' {dtAu{t - r)} • Au{t - r') ■ k['^\t - r') 
*^°° Jo Jo 

= lim [ dr [ dr' {dtA^^it - r)} ■ A,^{t - t') ■ K^^\t - r') = , (37) 

*^°° Jo Jo 

hence leading to the second expectation value vanishing. Along the same line, we can find, 

too, that the last expectation value, limt_j.oo Tr[{Pi(t), (5i(t)}+ ■ p(0)] =0. As a result, the 



steady-state heat current given in (|8pL|) reduces to 

K 






(iv) 



lim Tr 

2M t^oo 



Pi(t),g2(t) ■/)(o) 



(38) 



By applying the same technique, we can also arrive at the expression 



J. 



(N) 



K 



lim Tr 

2M t^oo 



Piv(t) , Qiv-i(t) ■p(0) 



(39) 



We can then find that limj^oo Qjv-i(^) appears directly from limt_>oo Q2{t) given in (34b ) with 



substitution of both A21 — ?■ Ajv-1,1 and A2M — ?■ Ajv-i,]v as well as \im.t_^oo -Pjv(^) comes from 



limt_>oo Pi{t) given in (34c) with An — )■ A^i and Aij^ — )■ Ajvat. Further, Eq. (30) straightfor- 



wardly gives rise to ^iv-i,i(s) = ^2iv('S) and ^jv-i,iv(s) = ^2i('S) as well as Ami{s) = ^ijv(s) 



and Ann{s) = ^11 (s). Substituting all these results into (39) and then with the aid of (36), 
we can verify that J^out = —J'in indeed. 

Now we are ready to derive an explicit expression of the steady-state heat current, based 
on the results obtained from the previous paragraphs. Then it turns out that (cf. Appendix 



B 






(,) _ hu'.Y. f^\'-' 



TT 



M 



2Af+l 

i=o 



K{-Zj) ■ hf,{zj) 



(40) 
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which is, in fact, vahd for arbitrary values of the input parameters {Q, K,,uJd,'yo)- Here the 
primed sum denoted by ^ . [■ ■ ■ ] means that if one of Zj^s is repeated, say, Z2m = ^2]v+i and so 
hj^{s) = (s + Z2js,)'^ ■ /]v(s) where /jv(s) = Ilj^o^ (s + Zj) and /i'^(— ;Z2jv) = 0, then it is needed 
to consider this primed sum spht into two parts such as X]i=o^ [' ' ' ] + A(2A^, 2N + 1), where 
h'j^{—Zj) = (— Zj + 2;2]v)^-/^(— -2j) 7^ 0; the extra part denoted by A(2A^,2A^+1) is contributed 
solely by the multiple root, s = —Z2n [cf- (41 3)-(41:)]. Fig. It] plots the typical behaviors 



of hj^{s). We also point out that the expression of heat current in (40), valid for A^ > 3 



r(2) 



corresponds to the heat current J'^l given in (25); note, however, that the denominator of 



each summand is here given by h'j^{—Zj) ■ hj^{zj) whereas it is in form of h[{—Zj) ■ hi{zj) for 

N = 2. 



To simplify the expression in (40), we substitute (18) into this and then apply the same 



technique as for A^ = 1, 2. Then we can straightforwardly obtain the exact closed expression 



(N) 



<-/m 



where 



'M2 



(i.) 



1 1 



2Af+l 



^ 2hwW. ^^yiV-2 ^, z]-{^,Xh)-^^Ah)') 



TT 



M 



j=0 



Ki-Zj) ■ h^{zj) 



(i^) 



Mio 



m) 



2]V-4 



2N+1 



2-^ h' 



Ki-Zj) ■ h^{zj) 



(41) 



m) 



note that this is in the unit of 1/ (frequency)^ as is the case for (I2) given in (27i). Here we 
also employed the sum rule given by ^ . z'j/{h'j^{—Zj) ■ /ijv(-Zj)} = for n odd (cf. Appendix 



r{N) 



E|). If Z2n = Z2N+1, then the heat current J7in contains the terms contributed solely by this 
licitly given by J7a 



multiple root, exphcitly given by J7a = J a (Pi) — Ja (Pn), where 



t(Jv) 

•J A 



(/3p 



Here /i = 1, A^, and /jv(— ^2? 



^^A\^ 



w„ 



27r ■ /iv( — ^2iv) ■ /iv(2;2jv) ■ Z2n 

= /i'^(-2;2jv)/2, as well as 
1 



X 



^(^) 



(/3. 






1 



(^d 



2N-1 



A^ 2 

i=o ^J 



"2jv 



^2jv 



[W/.(^2iV - (^d) - '^Z2NUJd{^Pf,{z2N) - V',9p(Wd)}] 



where u^ = 27r//3^/i, and the trigamma function ip^ {z2n) = '^^^\z2n / ^ ^i) ■ If one of Zj^s is 



lb) 



(41 3) 



repeated over than twice, we can straightforwardly generalize the result given in (41:) with 
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the help of (BIO) and (B14)-(B15) [29]. It is also worthwhile to mention that the compact 



expression given in (41 ) can be rewritten in terms of the input parameters (fi, k, Ud, 7o) with 



the aid of all coefficients of h^{s) explicitly given in (D3 ) and their symmetric properties like 



given in (11 ) used for A^ = 1, 2; the resultant expression will, however, be highly complicated 



even for A^ = 3. Eq. (41) is, in fact, the central result of this paper. 

Next we consider the semiclassical behavior of the heat current jTi^ . To do so, we expand 



the digamma function and its derivative given in (41)-(4lF) [cf. Appendix [B] , which will, in 
the semiclassical limit of /3^/i — )■ 0, gives rise to 






(iv) 



j}r\fi') + Ji^\h^) + JJ^^Kh^) + 0{h'). 



Here the leading term is given by the classical heat current 



(^)('hON 



Jtr{n 



fi: 7c 



(i») 



kP ' " Vft P^. 

If Z2n = ^2jv+i, then this classical value contains JTaVIi = JjCciWi) — ^i^i (/3]v), where 



(42) 



(43) 



r(iv) 



J^ >(B 



-co^^in/M) 



2N-2 



2N~1 



2/3^ ■ fN{ — Z2N) ■ fN{Z2N) ■ Z2p 



Wd - 6 {z2m - Wd) 1 + ^ 



"2iv 



z2 - 5-2 



m 



The first quantum correction is given by 

1 



J!:^\f^') 



/ K \ 2JV-2 
fc2 4 2 1'^ 



2N+1 



(/3i -MJ2 i;i 



Ki-^j) ■ hN{zj) 



(44) 



if necessary, with 



An) 



r '(B 



24/jv(-Z2iv) ■ fN{z2N) 



2N-1 



Z2i^ 



+ 2 {Z2r, - Ud) [l + J2 



^2n 



j=o ^i ^2^ 



W) 



The next quantum correction in Oif?^ is non-vanishing only if Z2n = Z2n+i in such a way 
that 



2Af--l 



^2n + \^2n 



U1.I 



i+E 



"2]V 



j=0 ^i ^2iv 



.(.).. . ^ 2 C(3)n3a;^ 72(^/^)2.-2 ^2 

' (45) 

where the symbol C^in) denotes the Riemann zeta function. In fact, all higher-order quantum 
corrections in closed form will straightforwardly come out. 
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It is now interesting to directly compare the classical result given in (43) with Fourier's 
law of heat conduction. This allows us to identify the classical heat conductivity Kp as 
(A^k^7q/M^) ■ (Ijv), which depends on chain length A^ hence violating Fourier's law already. 



From the same comparison of the quantum result given in (42)-(45), we can easily find 
the "effective" heat conductivity, which depends even on temperature due to the quantum- 
correction contributions. Therefore, we may argue that non-universal behaviors of the (effec- 
tive) heat conductivity, especially in low- dimensional lattices lying in the low-temperature 
regime (not only the harmonic chain under our investigation, as briefly stated in Sec. [l]), are 



ascribed by the non-classical contributions, as explicitly given in (44)- (45) for the harmonic 
chain, which are, in fact, not proportional to Ti — T^^ any longer. 

for 
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The behaviors of heat current jTJn versus chain length A^ are plotted in Figs. |8 
various input parameters. First, it turns out that in the weak-coupling regime imposed by 
k/M, 7^ <^ Q"^, the low-magnitude heat currents are typically acquired, as expected from 
the results for A^ = 1, 2. They also reveal the almost A^-independent behaviors (for A^ > 2). 
This can be understood from the forms of the matrix elements ^ii(s) and Ai2{s) in the 
weak-coupling limit; with the aid of sin A^0 = sin0 ■ cos(A^ — 1)0 + cos0 ■ sin(A^ — 1)0, we 



can exactly rewrite ^ii(s) in (30) as 



6 — c ■ cos — c ■ sin ■ cot(A^ — 1)0 

(46) 



^2 _ ^2 _|_ 2c ■ (c COS (f) — 6) ■ COS + 2c ■ (c COS — (5) ■ sin ■ cot(A^ — 1)0 ' 
where cot(A^ — 1)0 = —i coth(A^ — l)r, expressed in terms of the real number r = —i(f) > 
with 



coth(Ar - l)r = ^ _ (^!.)2(.-i) - 1 ■ &) 



In the weak-coupling limit leading to \c\ <^ b + 2c, Eqs. (32z) and (321) allow us to easily 
have 

which gives rise to coth(A^ — l)r — )■ 1 for A^ large enough (a fairly good approximation 
even for A^ = 3). Accordingly, the matrix element ^ii(s) reduces to be A^- independent. 
Along the same line, we can do the same job for Ai2{s),Ain{s) and A2n{s), respectively. 
Consequently, the heat current J7in reduces to be A^-independent in this regime. In this 
context, it is also worthwhile to mention that this behavior of heat current is consistent to 
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the result by Asadian et al. in [T], which was obtained from the consideration of a harmonic 



chain restricted to the rotating wave approximation of the isolated chain H^ given in (li) 
as well as to the Born- Mar kovian regime imposed by the weak-coupling and the Ohniic 
damping [uod -^ oo). Their result for steady-state heat current can be rewritten in terms 
of our notation, i.e., with V — ?■ fi;/2 as well as ri,rjv — > 7o in their equation (25), as the 
A^-independent expression 

,2 



J^H = ^ (Ib-m) ■ hn {{h),^ - (n),J , (47) 



where (Ib-m) '■= (1/2) {(^/^)^ + (^7o)^}~^; and the average excitation number (n)a^ = 
l/(e'''''*^ — 1). In the figures, this approximation is compared with our exact result denoted 
by J7in . It is then shown that this may be a good approximation in the weak-coupling 
regime, as expected, whereas it is not case beyond the weak-coupling regime. 

Next we pay attention to the above behavior of heat current beyond the weak-coupling 
regime, which has so far not been systematically explored. As demonstrated in the figures, 
the heat current increases with increase of the intra-coupling strength n for a given chain- 
bath coupling strength characterized by the imposed damping parameter 7^, and reaches its 
maximum value at some specific coupling strength Kr "resonant" to the chain-bath coupling 
strength. With further increase of the intra-coupling strength, the heat current decreases 
very slowly, whereas this behavior cannot be found from the Born-Markovian result given 



in (47). Also, the heat current typically behaves in such a way that its magnitude is at 
the maximum with A^ = 1, and then gradually decreases with increase of chain length A^, 
being in fact almost A^- independent in the range of A^ large enough. This may already be 



qualitatively underwood from the behavior of coth(A^ — l)r given in (46) with respect to A^. 
As a result, Fourier's law proves violated also in this regime. 

VII. CONCLUSION 

In summary, we rigorously derived an exact expression of the steady-state heat cur- 
rent through a chain of quantum Brownian oscillators coupled to two separate baths. It 
was obtained indeed for arbitrary coupling strengths both in the intra-coupling between 
two nearest-neighboring chain elements and in the chain-bath coupling, and in the Drude- 
UUersma damping model in order to look at the behavior of this heat current beyond the 
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Born- Mar kovian regime. Then we systematically observed that in the weak-coupling regime, 
the heat current with its low-magnitude simply reduces to be almost independent of the 
chain length. In the regime beyond the weak-coupling, on the other hand, we found the 
possibilities for the heat power amplification by observing the increase of the heat-current 
magnitude. 

By this study we could explore the relevance between the coupling strengths as input 
parameters and the magnitude of the output heat-current, which may be considered a fun- 
damental issue for building a quantum thermodynamic engine with high power. We also 
believe that our finding will provide a useful starting point for the analytical approach to the 
steady-state heat current through a more general type of quantum Brownian chains beyond 
the weak-coupling regime. 
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Appendix A: Derivation of correlation function in Eq. (3a) 



To derive a closed form of the correlation function in ([s]), we first employ the identity 



where the so-called Matsubara frequencies Un = 27m/ {(3h) |15| . Substituting (Al) into ([s]) 
and then applying the integral identities 



°° cos{ay) _ n e-"'' r°° dy cos{ay) _ it b e-"^ - c e-"'' 

^ y^ + b^ ~ 2~b~ ' io W+W(!lF+^ ~ 2 6c(62-c2) ' ^ ' 



we can straightforwardly obtain 



/ du — ^ 2 coth I — — I cos|a;(t — t )| 

Jo ^ +^d V 2 / 



TT 



e 



-uia l<-t 



(3hwd /3h ^ 

n=l 



' + |E "' ,._T ■ (A3) 



K-u} 
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We next apply to this expression both sum rules 

El 1 — Try cot(7ry) 

and 



n=l 



n^ — y^ 



2^/ 



.,,2 



2n p-«" 



-' n^ — y^ 



expressed in terms of the Lerch function [ 24j 

oo 



n=0 



n + f)' 



which finally gives rise to the closed expression in ( 3 1 ) 



(A4) 



(A5) 



(A6) 



Appendix B: Derivation of Eqs. (17)-(20a) 



We first substitute (13) into (16), which immediately yields 



.{d), 



A(t) , i,At)]^ = 2M/hn ^ dr |^Xd(t - r) \ K^it - r) 
Here the bath correlation function is explicitly given by 

,2 



(Bl) 



Klf\t 



n = ^^^4^Kcot(M^]e--l- 



27r 



*-" + E 



2^.g-na;M|t-r| 



n=0 



{n + Ud/uj^) ■ {n - UJdli^i}i 



[cf. (|A6|)]. We can easily evaluate the integral in (|B1|) explicitly, which leads to 



Pi{t)AM 



W huj^^M 



TT {Zo - Zi){zi - Z2){Z2 - Zo) 



Here 

y^ = ^ 

where 



n^) = E 



^i+i ^i+2/ ^i 



u + Zi 



By means of the identity of the digamma function [;25j 

A^ (n J- n\ ( 



n=0 



{n + a) {n + b) a — b 



(B2) 



(B3) 



cot (^) .r(c.,) + f: (—V + ^V) ■^(^"^) ' (^4) 



(B5) 



(B6) 
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we can easily rewrite the sumniation in (B4) as 

2 






2tt ^ I 27r 



[Ud -^ -UJd] 



(B7) 



Substituting now the expression in (B4) into (B3) and subsequently into (16), we can finally 



arrive at the result given in (17). 



Next let us derive the expression of heat current given in (20) in terms of the input 



parameters {Q,uJd,'yo) only, expanded in the semiclassical limit. To do so, we first plug into 
^ the expansions given by cot{y) = ^^^^ (-1)" {22"/(2r2)!} 52„y2""^ for < |y| < vr 
and ipiv) = ~^/y ~ 7e + X]^i (~1)"^^C(^ + l)l/"j here the Bernoulli numbers i?2n; the 
Euler constant 7e = 0.5772 ■ ■ ■ , and the Riemann zeta function ({n + 1), as well as i?2n = 
2 (— l)"'^^{(2?7,)!/(27r)^"} C(2n) [24j. After some steps of algebraic manipulation, this gives 
rise to 



X:' = Jl"ih') + J.l^Kh') + J!iKh') + Jl^lKh') + J^^lKh'') + J!ilKh') + om . (bs) 



Here we have the leading term 



/3i /3i' 



j=0 



Zj [Zj 



UJd) 



i^j ~ ^i+i)(^i ~ ^i+2)(^j + '^d) 



(B9) 



To evaluate this summation explicitly, we take into account the technique of partial fraction 
for two polynomials P{s) and Q{s) with degQ{s) < deg P(s) = n, where P{s) = {s — bi){s — 
62) ■ ■ ■ (s — bn) with bj 7^ bk for j 7^ k. This is explicitly given by [22] 

Qis) ^ Q{b,) 



m :-- 



Pis) § P'ib,) ■ is - b,) • 



(BIO) 



Applying this relation, the summation in (B9) easily reduces to 



S {S + Ud) 



(Bll) 

[S + Zo) [S + Zi) [S + Z2' 

which allows us to have the classical result in (20i]). Here we also used the relations in (11 ) 
Next, we consider the quantum corrections 



J'^Kn-) 



J}Pi^') 



^^-^^"(A-AOEr^ ^K^.--) 
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~^ i^i ~ ^i+i) i^i - ^i+2) i^i + ^d) 



25 . 32 ■ 5 



{PI-PDUtt^ 



^] (4 ~ ^rf) 



j=0 



(^i " ^i+i)(^i ~ ^i+2)(^j + '-^d) 



(B12) 
(B13) 
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Here we used ({2) = tt'^/Q and ^(4) = vr^/QO. Applying again (BIO) to these two summations 



and then evaluating them at s = Ud, respectively, we can finally arrive at the result in (20) 



In fact, every quantum correction with the odd-degree h-powei in (B8) is shown to vanish 



indeed by applying the same technique. Finally we point out that if one of the roots b^, is 



repeated m times in (BIO), then the expansion for f(s) contains the terms of form 

Ai 



s - 6. 



+ 



where 



lim 



A2 

is - K)^ 



1 / d 
r! V ds 



+ 



+ 



Xr 



{{s-Kr-fis)} 



This will be used in Sec. IVII 



(B14) 



(B15) 



Appendix C: Evaluation of Eqs. (23)-(25) 



We consider the double integral given in (23) 



(I2) := hm fdr f dr' Xiit - r) ■ dt{x2{t - r')} ■ K\^\r - r') . (CI) 

*-^°° Jo Jo 

We substitute ( 12^ ) and (B2) into this. In doing so, let Xd{t) -^ Xi(^) expressed in terms of 

^ij's, and Xd{t) -^ X2{t) expressed in terms of Z2./s. Then it turns out that 

hcol % Z^l,k=0 (^l,i+l " ^l.j+a) (^2,fc+l ~ ^2,^+2) ^2,A: 



(I. 



TV cot 



where / = 0, 1, 2, and 



27rM {(2:1,0 - ^1,1) (^1,1 - ^1,2) (2:1,2 - 2:1,0)} ■ {{zi^i -^ Z2,i)} 
' I3^f\hjd\ .- , -^ 2n 



X 



(I 



2,lJi^d 



+ 5: 

n=0 



{n + Ud/ujp,) ■ {n - ujd/ujf,) 



(I 



2,l)nuiu ( ) 



(C2) 



*-^°° Jo Jo 



1 1 

+ 



(C3) 

Zl,j + Z2,k \ZlJ + a Z2,k + CiJ 

Here we also used the integral identity given in (15). 
Next we consider 

(II2) := lim f dr [ dr' d^ixiit - r)} ■ X2{t - r') ■ Kf{T - r') . (C4) 

*-^°° Jo Jo 

Along the same line, this reduces to the expression given in (C2) but with exchange of Zij 



and Z2,k- Noting Xi(0) = X2(0) = from (12i), we can first find that (I2) + (II2) vanishes 
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indeed, and so does Eq. (23). Next, (I2) — (II2) gives rise to an explicit evaluation of the 



integral in (24) and subsequently the exact result in (25) expressed in terms of {zj\ Zj = Zij} 



and {Zj\zj = Z2j}. 



Appendix D: Mathematical supplements for Eq. (33) 



First, let us acquire an explicit expression of den(s; A^) which leads to the polynomial 
hN{s) = den(s; A^) ■ [s + UdY- To do so, we rewrite Eq. (32") as F,^{s) = Fi{s) ■T,v(s), where 

N~l 00 



.2\J 



T.is) = J2^^(')y"'-ws)r = Yl (27 + 1) ■^"^''"'- (^'-^"') 



u=0 



3=0 



N 



(Dl) 



With the aid of (32a)- (32b), we can then find that 



Aii{s) 
AAs) 



2^- [T,{s) + 2{s7(s) + c}-T,_i(s)]/den(s;iV) 



-,2-JV 



(-c) [T^-i{s) + 2{s7(s) + c}-T^_2(s)]/den(s;Ar) 



Ai^{s) = (-c)^- Vden(s; AT) , ^^(s) = a (-c)^-Vden(s; A^) . 



(D2a) 
(D2b) 
(D2c) 



Substituting (Dl) into (32b) and then applying Pascal's rule (^) = ( ^ ) + (fc_i)) we can 



finally arrive, after some algebraic manipulations, at the expression 

den(s;A^) = 2 
'N -I 



-,1 — N 



£.r"r'){(r.:)[''^-<-'-»''-^^]- 



[(f + 2{x + (-c)} d + A (-c) x] 



X 2" 



k+2n jN—k—n—2 



{-cY+'\ (D3) 



where x{s) := 37^(3) and d{s) := s^ + fi^. From this, all coefficients of den(s; A^) can exactly 
be determined, which are non-negative, as shown; e.g., the highest s-power term is given by 
g2N+2^^g _l_ ^^^2 ^^^ ^]^g second highest s-power term is 2ci;rfS^^+Y(s + w^)^- Substituting 



(D3) into (D2a)-(D2c), we can get the explicit expressions in (33). 



Next, we let us prove that R,e{zj) > for all z/s satisfying h 



N\ ^3, 



0, which is needed 



for applying the final value theorem of the Laplace transform. Due to the nonnegativeness of 



all coefficients given in (D3), it suffices to prove that den(s; A^) 7^ for any purely imaginary 
number s = ir, where r G M. We assume that den(zr; A^) = 0, though. Then, its conjugate 
number s = —ir should also satisfy the equality, den(— ir; A^) = 0. By applying these two 
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equality conditions to (30o) simultaneously, we can obtain both 

c ■ sin N(f) = 



r^ + u;^ 



0+"^^^^] ■sin(iV-l)0 



c c + 






■sin(A^ + l)0 = 5 <^2 c + 






5} -sin N (p. 



(D4a) 



(D4b) 



From this, we notice that sin A^0 7^ if sin0 7^ 0. Combining (D4a) and (D4b) to eliminate 
the sine functions therein, we can easily acquire 



.2 , ,4 



{{^,Ud + cy + c'}X' + ui{{^,Ud + cy + 3c'}X + 2c'u 







(D5) 



where X := r^ > 0. Then we see that each root Xg of this quadratic equation would have 
to meet the condition Ke{Xg) < 0, which, however, contradicts itself. Consequently, we 
cannot have any numbers Zj's being purely imaginary. It is, however, nontrivial indeed to 
extract all individual roots {—Zj)'s of the polynomial /liv(s) expressed explicitly in terms 
of the parameters {Q, K,Ud,'yo), even for A^ = 3 when we need to deal with the 8th-degree 
polynomial h^ls). 



Appendix E: Evaluation of Eqs. (35)-(40) 



To evaluate the integral in ( 35 ) explicitly, we first consider the double integral 

ft 



(U, := lim / drf{t-r) [ dr'g{t-r') 
*^°° Jo Jo 



■ e 



-a \t—t'\ 



(El) 



where f{t) = An{t) and g{t) = A2i{t). By applying the technique in (C3) used for A^ = 2 



we can transform (El) into 



(U)„ = ^[l dr[fit)-gir) + g 



W-/(r)]e"^j(a). (E2) 

Let /(s) = £{f{t)}{s) and g{s) = C{g{t)}{s). Now we consider the product rule, which 
reads as £{/i(t) f2(t)}{s) = (l/27ri) f^\-^ du fi{u) f2{s — u) j22]; here the integration is car- 
ried out along the vertical line, Re(M) = Ci that lies entirely within the region of convergence 



of fi{u). Then we can easily rewrite (E2) as 

"^1+^°° du 



(I. 



C\—lOO 



u — a 



{f^u)-g{-u) + -g{u)f{-u)} , 



(E3) 



where f{u) = uAn{u) and g{u) = A2i{u). Here we also used >C{/g dr g{T)}{s) = g{s)/s. 
With the aid of ( 32|i ), the integrand given by {u [Aii{u) A2ii—u) — A2iiu) Au{—u)]} can 
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be transformed into {—u[Ain{u)A2n{~u) — A2n{u) AiNi—u)]}, which immediately allows 



us to obtain the relation given in (36) 



Therefore we can evaluate the integral in (El) by plugging f{t) — )■ Ai,^{t) and g{t) — )■ 
^2iv(^) giving rise to — (Ijv)a; in fact, Ai,^{t) and A2N{t) are simpler in form than Aii{t) and 
A2i{t), respectively. First we rewrite the expressions in (p3| as 



2Af+l 



Aim{s) 



L^ hi 



-Zj + UJd) 



2N+1 



h'(—Zj. 



is + Za 



A (A - ( A^-^ \^' (~^J +^d) ■ hi{-Zj] 



(E4) 



respectively. The meaning of the primed sum denoted by ^ . is explicitly given below (40) 



which must be treated with care if one of z/s is repeated [cf. (BIO) and (B14)-(B15)]. Then 
it easily follows that 



fit) 



git) 



2N+1 

1"-' E' 

j=0 

2N+1 

j=Q 



Ki-Zj) 

'-Zj +ujd) ■ hii-z 



Ki-Zj) 



1) ^-Zjt 



(E5a) 
(E5b) 



each of which, in the Zfc-degenerate case, contains the terms resulting from (B14)-(B15) in 
such a way that 



£-^ 



'^m—r 



is + Zk)' 



it) 



Xr 



f 



m—r—1 



ira — r — 1)\ 



-Zkt 



(E6) 



where r = 0, 1,--- ,m — 1. We now substitute this result into (El ) and evaluate the double 



integral explicitly, which is the same in form as the integral in (C3) considered for chain 
length N = 2 only. Therefore we can straightforwardly obtain 



2Af+l 



(I. 



-c 



,2n-3 



2-^ h' 



I Zj ■ {zj - Ud)^ ■ hi{-Zk) 



j,k 



Ki-Zi) ■ Ki-Zk) ■ izi + Zk) \zj + a 



1 



1 



Zk + a 



(E7) 



Applying the technique of partial fraction given in (|B10|), this can be simplified as 

2N+1 



4 ■ (4 - ^i) 



(E8) 



-ivv -Zj)-Kizj)-ia + Zj) 
This allows us to have an explicit evaluation of the integral in ( p5| and then that of the heat 
current 



7^ 

<-/in 



(N) 



cot 



hhjA 



2k^T^ 
n ■ (Iiv)„ 



cot 



hWri 



Z rVTll Aj 



ih 



'^d 



+ 



2M 

oo 

^^ {in + U!d/uJi){n-UJd/uji) (n + ujd/ujN)in - Ud/uj^) 



Ul 



n ■ (In) 



nujff 



(E9) 
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as provided in ( 40 ) . 



Next let us prove the sum rule given by 

2Ar+l 



V' -^-, x^ , . X = (ElO) 



Ki~Zj) ■ h^{zj) 



for n odd, which is used for (41 ). First we rewrite hj^(zj) as nA;=o (^j "*" ^k) = Y[k=o i~^j + 
Z2N+2+k), where we introduce Z2M+2+k '■= ~Zk. Then it turns out that h'^{—Zj) ■ h^{zj) = 
Ylt1t^{-Zj + Zk) where k ^ j. Next let ifjv(s) := h^i^s) ■ h^i-s) = YllZo^{s + Zk), and we 
consider 

^^'^ '" ^M^ ^ £ HU-Zk) ■ {s + Zk) ■ ^^^^^ 

Then we can easily obtain 

2Af+l „ 

F(0) = = (-1)"+^ {1 - (-1)"} 5^' „,^l, . , , (E12) 



j=0 ^'n(-^j) ■ ^n{Zj) ' 



which immediately gives rise to the sum rule in (ElO). In case that one of 2;j's is repeated, 
it is also straightforward to verify this result. 
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X 



FIG. 1: 



Fig. ll (Color online) y = /C}j (bath correlation function) versus x = t (time), given in 

1; solid line plotted at T^ = 1 (low temperature) 



(3i). Here we set h = k^ 



M 



% 



while dashed line at T^ = 10 (high temperature). From top to bottom at x = 0.75, 1st: 
(green dash: Ud = 1); 2nd: (blue solid: Ud = 1)', 3rd: (red dash: Ud = 10); 4th: (black 
solid: Ud = 10). For the 1st, 2nd and 3rd lines, we have w^ < w^ = 27rT^. For the 4th, 
on the other hand, we have Ud > u)^, which gives rise to the low-temperature behavior of 
JCjj, (t) characterized by appearance of its negative- valued region; all four lines diverge at 
t = anyhow due to the result given by JCfj, (0) = kj^T^M'^^tUdil + 2aX]^i V(^ + '^)} 



with a = hud/{2TT ksT^), directly obtained from (A3). 
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FIG. 2: 



Fig. 2 (Color online) y = /C}j (bath correlation function) versus x = T^ (temperature). 
Here we set /i = /cb = M = 7^ = 1; solid line representing the typical early-time behavior, 
plotted at time t = 0.3 while dashed line representing the late-time behavior, plotted at 
t = 2.5. From top to bottom at x = 5, 1st: (green solid: Ud = 1); 2nd: (blue solid: 
Ud = 10); 3rd: (red dash: Ud = 1); 4th: (black dash: Ud = 10). As demonstrated, the 
correlation function has no singularities at Ud/oo^ = 1, 2, ■ ■ ■ , where w^ = 2'kx. 
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FIG. 3: 

Fig. 3: (Color online) y = jTjJ; (heat current) versus x = Ti (hot-bath temperature), in 
the low-temperature regime where the cold-bath temperature is imposed by Ti/ = 0.1. As 
such, we see that y = at the thermal equilibrium point, x = 0.1. Here we set h = kj^ = 
M = Q = 1, and Ud = 10; solid line plotted for the quantum- mechanical heat current given 
in ((l9|) while dashed line for its classical counterpart in (20 1). From top to bottom at a; = 1, 



1st: (black dash: 7^ = 1); 2nd: (green solid: 7^ = 1); 3rd: (blue dash: 7^ = 0.2); 4th: (red 
solid: 7o = 0.2); 'j^ = 0.2 represents the weak coupling 7^ ^ fi between the single oscillator 
and two baths. 
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FIG. 4: 



Fig. |4J (Color online) The same plot as in Fig. [sl in the high-temperature regime where 



the cold-bath temperature is imposed by Ti/ = 2. As such, we see that y = J-, 



(1) 



at 



the thermal equilibrium point, x = 2. In the region of x > 2, the quantum-mechanical 
and its classical values almost overlap each other; on the other hand, y < for x < 2. 
From top to bottom at x = 0, 1st: (red solid: 7o = 0.2); 2nd: (blue dash: 7^ = 0.2); 3rd: 
(green solid: 7^ = 1); 4th: (black dash: 7^ = 1). It can further be verified that in case 
that the cold-bath temperature is even higher, i.e., given by Ty > 2 ("classical regime"), the 
quantum-mechanical heat current will more strongly overlap its classical counterpart. 



34 



0.15- 



0.10- 



0.05- 




X 



FIG. 5: 



r(2) 



Fig. 5 (Color online) y = J7in (heat current) versus x = Ti (hot-bath temperature) 



given in (26), in the low-temperature regime where the cold-bath temperature is imposed 



by T2 = 0.1. As such, we see that y = at the thermal equilibrium point, a; = 0.1. Here we 
set h = kj^ = M = Q = 1, and Ud = 10. Solid lines plotted for 7^ = 1, from top to bottom 
at X = 1.5, 1st: (green: k = 1); 2nd: (khaki: k, = 0.5); 3rd: (red: k, = 0.2). Dashed lines 
for 7o = 0.2, from top to bottom at a; = 1.5, 1st: (blue: k, = 1); 2nd: (brown: k = 1.5); 3rd: 
(black: k, = 0.2). The 3rd dashed line represents the weak-coupling regime, k/M, 7^ ^ Q^. 
Notably, this is of higher value than the 3rd solid. Also, for 7^ = 0.2 the maximum (or 
resonant) heat current is obtained at k, = 1 = k,^ while with further increase of /t > k,^, the 
current decreases very slowly. 
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FIG. 6: 
Fig. pj (Color online) The same plot as in Fig. M in the high-temperature regime where 

(2) 

the cold-bath temperature is imposed by T2 = 2. As such, y = jTJn = at the thermal 
equilibrium point, x = 2. In the region of a; > 2 the heat current reveals the behavior of its 
classical counterpart, being proportional to Ti — T2, while for x < 2 it is negative- valued. 
As shown, we have the same lines for x > T2 as in Fig. |5| except that the 1st top dashed 
line (at x = 10): (brown: 7^ = 0.2, and k, = 2.2 instead of 1.5), and 2nd top dash: (blue: 
7o = 0.2 and k, = 1), i.e., for 7^ = 0.2 the maximum heat current appears at k = 1 = Kr 
while with further increase of k > k^, the current decreases very slowly. 
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FIG. 7: 



Fig.uj (Color online) y = hpf{s)- 10 ^^ versus x = s, given in (33). Here we set ^2 = M = 1 
and Ud = 10 as well as A^ = 5. From top to bottom at x = —7.5, 1st) h^^i{s) ■ 10^^"*^: (red 
solid: 7„ = 0.2 and k = 0.2) with a single multiple root, s = -9.79840398; 2nd) /i5,2(s)-10-": 
(blue dash: 7„ = 0.2 and k = 1) with a multiple root, s = -9.80006223; 3rd) /i5,3(s) ■ 10"": 
(green solid: 7„ = 1 and k = 0.2) with a multiple root, s = -8.89222703; 4th) h^^is) ■ 10"": 
(black dash: 7^ = 1 and k = 1) with a multiple root, s = —8.90443052. All the multiple 
roots have degeneracy m = 2. We also have /i5,i(0) = 0.00000361 ^ 0; /i5,2(0) = 0.00005500; 
/i53(0) = 0.00000361; /i5,i(0) = 0.00005500; in fact, /i^(0) increases with increase of A^. 
The numerical analysis of the exact expression of hj^{s) reveals that there is no additional 
complex- valued multiple root of the above four functions. These properties of hpf{s) are 
verified to be valid for many different choices of (7^, k), and A^ = 3,4, ■ ■ ■ ,20. 
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FIG. 8: 

(Color online) y = J7in (heat current) versus x = N (chain length), given in 
((41)-( 4lp) )), in the low-temperature regime imposed by Ti = 1.1 and Tjv = 0.1. Here we 
set h = kj^ = M = Q = 1, and Ua = 10. Dashed lines, from top to bottom at x = 20, 1st: 
(red with solid circles: 7^ = k = 1); 2nd: (green with diamonds: 7^ = 0.2 and k = 1); 3rd: 
(black with diamonds: 7^ = k = 0.2); 4th: (blue with solid circles: 7^ = 1 and k = 0.2). 



Fig. 



r(iv) 



In comparison, solid lines plotted for J'b-m in given (47), from top to bottom, in the same 



order as for the dashed lines. As demonstrated, this represents a good approximation in the 
weak-coupling regime. 
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FIG. 9: 



Fig. PI (Color online) The same plot as in Fig.^ in the high-temperature regime imposed 
by Ti = 3 and Tjv = 2. 
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FIG. 10: 



Fig. 10 (Color online) y = jTJn (heat current) versus x = k (intra-coupling strength), 
given in ((41)-(4l3)). Here we set h = k^ = M = Q = 1, and Ud = 10, as well as 7^ = 0.2. 



Dashed lines, from top to bottom at x = 1, 1st: (red with solid circles: A^ = 3 as well as 
Tpf = 3 and Tjv = 2, with its maximum at x = k,^ = 2.2); 2nd: (blue with diamonds: A^ = 6 
as well as Tjv = 3 and Tjv = 2, with Kj^ = 2.8); 3rd: (green with solid circles: A^ = 3 as 
well as Tjv = 1.1 and T^ = 0.1, with Kr = 1); 4th: (black with diamonds: A^ = 6 as well as 
T,^ = 1.1 and T,^ = 0.1, with Kr = 1.4). This shape of the heat current with respect to the 
intra-coupling strength is verified to be true for different choices of the temperature range 
and all other input parameters. In comparison, two solid lines plotted for JTb.m in given (47): 



the red upper for Tjv = 3 and T^ = 2, asymptotically approaching 0.09862324 with x — )■ 00, 
and the green lower for T^ = 1.1 and Tpf = 0.1, approaching 0.06746888 with x — t- 00. 
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